Utilisation des fonctions MH

Exemple simple pour retrouver la variances et la moyenne d’un échantillion gaussien

On utilise ici la version grand dimension et la version grand dimension parallélisé de l’algorithme Metropolis Hastings.

Cet exemple sert à montrer l’utilisation des fonctions MH.

note

Pour la version parallélisé il est demandé de fournir un écart type sd (possiblement un vecteur d’écart type) mais plus une fonction de transition.

n <- 10000# nombre d'observation
x <- c(3,6) #moyenne et variance
obs <- rnorm(n, x[1], sqrt(x[2])) #tirage dans une loi normale

transi = function(x) rnorm(length(x), x, 0.1) #fonction de transition entre x et xnew

loglik <- function(x, obs) sum( -1/2*log(2*pi*x[2]) -1/(2*x[2])* (obs-x[1])^2 ) #vraisemblence

#Metropolis hasting
x0 <- matrix(c(1,1), nrow = 1)
res <- MH_High_Dim(500, x0, transi, loglik, obs, verbatim = T)
plot(res, nrow = 2)

res_para <- MH_High_Dim_para(500, x0, sd = 0.1 , loglik, obs, verbatim = T, cores = 1)
plot(res_para, nrow = 2)

res %>% as.numeric()
## [1] 3.021748 6.034551
res_para %>% as.numeric()
## [1] 3.013567 5.900558

Note

On peut relancer l’algorithme tout en gardant les informations déjà compiler

res_para2 <- MH_High_Dim_para(100, res_para, sd = 0.1 , loglik, obs, verbatim = T, cores = 1)
plot(res_para2, nrow = 2)

# Exemple d'optimalité entre les 3 fct MH, MH_high_dim et MH_high_dim_para sur le modèle non lineaire mixte
m <- function(t, eta, phi) (phi[1] + eta)/(1+exp((phi[2]-t)/phi[3]))

#=======================================#
param <- list(rho2 = 0.1,
              sigma2 = 0.05,
              mu = c(5,4,0.5),
              omega2 = c(0.5,0.1,0.01))
              # mu = 5,
              # omega2 = 0.5)
F. <- length(param$omega2) #dimension de phi
#=======================================#

t <- seq(2,6, length.out = 10) #valeur des temps

G <- 40 ; ng = 4
dt <- NLME_obs(G, ng, t, param, m)
Y <- dt$obs

dt %>% ggplot(aes(t, obs, col = gen, group = id)) +
  geom_point() + geom_line() +
  theme(legend.position = 'null')

#=============================================================================#

Le but ici est de montre la différence de temps de calcul, a noté que certe MH est rapide ne fait pas la même chose que les deux autres fct.

De plus on montre

Phi <-  c(- 1/(2*param$sigma2), -1/(2*param$omega2), param$mu/param$omega2)

S <- function(phi)
  c(sum((Y - get_obs(m, dt, phi = phi) )^2 ),
    apply(phi^2, 2, sum),
    apply(phi  , 2, sum) )

#=============================================================================#
sd <- 0.02

transi <- function(x) rnorm(length(x), x, sd) #Fonction de transition pour MH_high_dim
transi.noopti <- function(x) rnorm(length(x), x, sd) %>% matrix(nrow = G) #Fonction de transition pour MH

loglik <- function(phi)
{
  stat <- S(phi)
  sum(Phi*stat)
}

phi0 <-  matrix(rnorm(F.*G, c(6,3,1), 0.1), ncol = G) %>% t
n.iter <- 1000

tic()
phi.MH <- MH(n.iter, phi0, transi.noopti, loglik, verbatim = T)
toc()
## 26.79 sec elapsed
tic()
phi.MH2 <- MH_High_Dim(n.iter, phi0 , transi, loglik, verbatim = T)
toc()
## 1063.157 sec elapsed
tic()
phi.MH3 <- MH_High_Dim_para(n.iter, phi0, sd, loglik, verbatim = T, cores = 1)
toc()
## 1057.408 sec elapsed
tic()
phi.MH4 <- MH_High_Dim_para(n.iter, phi0, sd, loglik, verbatim = T, cores = parallel::detectCores() -1)
toc()
## 46.532 sec elapsed
tic()
phi.MH5 <- MH_High_Dim_para_future(n.iter, phi0, sd, loglik, verbatim = T, cores = 1)
toc()
## 1295.729 sec elapsed
tic()
phi.MH6 <- MH_High_Dim_para_future(n.iter, phi0, sd, loglik, verbatim = T, cores = future::availableCores() -1)
toc()
## 91.745 sec elapsed
#=============================================================================#
#Affichage resultat
v <- attr(phi.MH, 'value') %>% lapply(function(x) as.numeric(x)) %>% as.data.frame %>% t %>% as.data.frame()

v %>% mutate(i = 1:nrow(.)) %>% melt(id = 'i')  %>% mutate(phi = rep(factor(c(1:F.)), each = nrow(v)*G)) %>%
  ggplot(aes(i, value, col = phi, group = variable)) +
  geom_line()

plot(phi.MH2, nrow = 2)

plot(phi.MH3, nrow = 2)

plot(phi.MH4, nrow = 2)

plot(phi.MH5, nrow = 2)

plot(phi.MH6, nrow = 2)

#=============================================================================#
n <- 800
m2 <- MH_burnin(phi.MH2, n) %>% attr('value') %>% melt(id = c('id','iter')) %>% group_by(interaction(variable,id)) %>%
  summarise(mean = mean(value), sd = sd(value))

m3 <- MH_burnin(phi.MH3, n) %>% attr('value') %>% melt(id = c('id','iter')) %>% group_by(interaction(variable,id)) %>%
  summarise(mean = mean(value), sd = sd(value))

m4 <- MH_burnin(phi.MH4, n) %>% attr('value') %>% melt(id = c('id','iter')) %>% group_by(interaction(variable,id)) %>%
  summarise(mean = mean(value), sd = sd(value))

m5 <- MH_burnin(phi.MH5, n) %>% attr('value') %>% melt(id = c('id','iter')) %>% group_by(interaction(variable,id)) %>%
  summarise(mean = mean(value), sd = sd(value))

m6 <- MH_burnin(phi.MH6, n) %>% attr('value') %>% melt(id = c('id','iter')) %>% group_by(interaction(variable,id)) %>%
  summarise(mean = mean(value), sd = sd(value))

mean(abs(m2$mean - m3$mean))
## [1] 0.5902355
mean(abs(m3$mean - m4$mean))
## [1] 0.07936126
mean(abs(m4$mean - m5$mean))
## [1] 0.07876204
mean(abs(m5$mean - m6$mean))
## [1] 0.0958769
mean(abs(m6$mean - m2$mean))
## [1] 0.5659555
mean(abs(m2$sd - m3$sd))
## [1] 0.02676858
mean(abs(m3$sd - m4$sd))
## [1] 0.01533202
mean(abs(m4$sd - m5$sd))
## [1] 0.01184366
mean(abs(m5$sd - m6$sd))
## [1] 0.01640094
mean(abs(m6$sd - m2$sd))
## [1] 0.0326134